function [fitCurve, fitStats] = cosFit(phase,mean,stdev,nPoints,dataFileName,sTitle)
 
% Set up figure to receive datasets and fitted curve
% =========================================================================
plotFig = clf;   % Clear current figure window
plotAxis = axes; % Create an axis object
figure(plotFig); % Create a figure graphics object


% Set some axis properties
% =========================================================================
set(plotAxis,'Units','normalized','OuterPosition',[0 0 1 1],'Fontsize',10);
set(plotAxis,'Box','on');
set(plotAxis,'XMinorTick','on');
set(plotAxis,'YMinorTick','on');
xlim = [0 1];
axes(plotAxis);
hold on;
 
% Plot the folded data points with stdev error bars
% =========================================================================
phase = phase(:);
mean  = mean(:);
stdev = stdev(:);
nPoints = nPoints(:);

plotHandle = errorbar(phase, mean, stdev, 'Parent', plotAxis,...
    'LineStyle','none', 'LineWidth', 0.5,...
    'Marker','o', 'MarkerSize',4, 'MarkerFaceColor','b');

%plotHandle = plot(phase, mean, 'Parent', plotAxis,...
%    'LineStyle','none', 'LineWidth', 0.5,...
%    'Marker','o', 'MarkerSize',4, 'MarkerFaceColor','b');

% Create a fit to the folded data points
% =========================================================================
fitFunc = fittype('A + B*cos(2*pi*x + phi)'); % equation to fit (A, B, phi)

fitOpts = fitoptions('method','NonlinearLeastSquares');
set(fitOpts,'Lower', [0.0 0.0 -pi])         % set the lower bounds on the coeffients
set(fitOpts,'Upper', [Inf Inf pi])          % set the lower bounds on the coeffients
set(fitOpts,'Startpoint', [0.2 0.5 0.0]);   % set the starting points
%set(fitOpts,'Weight', nPoints);            % weight given to point

[fitCurve, fitStats] = fit(phase, mean, fitFunc, fitOpts);

coeff = coeffvalues(fitCurve);
%err = confint(fitCurve, 0.99); % 3sigma
%err = confint(fitCurve, 0.95); % 2sigma
err = confint(fitCurve, 0.68); % 1sigma
disp('-------------------------------------------------');
disp(['A + B * cos(2 * pi * x + phi)']);
boundA_1 = coeff(1)-err(1);
boundA_2 = err(2) - coeff(1);
disp(['A = ', num2str(coeff(1)), '  +/- ', num2str(boundA_1), ' (',num2str(boundA_2),')']);
boundB_1 = coeff(2)-err(3);
boundB_2 = err(4) - coeff(2);
disp(['B = ', num2str(coeff(2)), '  +/- ', num2str(boundB_1), ' (',num2str(boundB_2),')']);
boundPhi_1 = (coeff(3)-err(5))/(2*pi);
boundPhi_2 = (err(6) - coeff(3))/(2*pi);
disp(['Phi = ', num2str(coeff(3)/(2*pi)), '  +/- ', num2str(boundPhi_1), ' (',num2str(boundPhi_2),') (x pi)']);
disp(['rsquared = ', num2str(fitStats.rsquare)]);

str1(1) = {'{\bfy = a + b cos(2 \pi t + \phi)}'};
str1(2) = {['a = ', num2str(coeff(1),'%2.4f'), ' (\pm ', num2str(boundA_1,'%2.4f'),')']};
str1(3) = {['b = ', num2str(coeff(2),'%2.4f'), ' (\pm ', num2str(boundB_1,'%2.4f'),')']};
str1(4) = {['\phi = ', num2str(coeff(3)/(2*pi),'%2.4f'), ' (\pm ', num2str(boundPhi_1,'%2.4f'),')',' \times 2\pi']};
str1(5) = {['r^{2} = ', num2str(fitStats.rsquare,'%1.3f')]};


fid = fopen('log.txt', 'a');
%disp([dataFileName]);
name=regexp(dataFileName,'p(\d*.\d*|\d*)-','tokens');
period=eval(cell2mat(name{1}));
phase=coeff(3)/pi;
phase_err = (coeff(3)-err(5))/pi;
a = coeff(1);
a_err = coeff(1)-err(1);
b = coeff(2);
b_err = coeff(2)-err(3);
r2 = fitStats.rsquare;
newLine = 'n';
msg = '';
msg = strcat(msg,sprintf('%.3f ', period));
msg = strcat(msg,sprintf('%10.4f %12.5f %12.5f ', a, a+a_err, a-a_err));
msg = strcat(msg,sprintf('%10.4f %12.5f %12.5f ', b, b+b_err, b-b_err));
msg = strcat(msg,sprintf('%10.4f %12.5f %12.5f ', phase, phase+phase_err, phase-phase_err));
msg = strcat(msg,sprintf('%10.4f', r2));
disp([msg]);
fwrite(fid, msg);
fprintf(fid, '\n');
fclose(fid);


% ========= Plot fit ==================
plotHandle = plot(fitCurve, 'fit', 0.95);
%plotHandle = plot(fitCurve,'predobs',0.95);
plotHandle = plot(fitCurve,'predfunc',0.95);
%plotHandle = plot(fitCurve,'residuals',0.68);
legend off;  % turn off legend from plot method call
set(plotHandle(1),'Color','black',...
     'LineStyle','-', 'LineWidth', 1.2,...
     'Marker','none', 'MarkerSize',6);
plotHandle = line([0.0 1.0], [coeff(1) coeff(1)]);
set(plotHandle,'Color','black',...
     'LineStyle',':', 'LineWidth', 0.1,...
     'Marker','none', 'MarkerSize',6);


%xlabel(plotAxis, 'Orbital phase /2\pi','FontName','Times New Roman', 'fontsize',12);
ylabel(plotAxis, 'mean speed /c','FontName','Times New Roman','fontsize',12);
%ylabel(plotAxis, 'Amplitude','FontName','Times New Roman','fontsize',12);
xlabel(plotAxis, 'Phase /2\pi','FontName','Times New Roman','fontsize',12);
title(plotAxis, sTitle,'fontsize',18,'FontName','Times New Roman');
text(0.75, 0.275, str1, 'HorizontalAlignment', 'left', 'fontsize',8);
%text(0.72, -0.65, str1, 'HorizontalAlignment', 'left');
ylim([0.23 0.28])
%ylim([-1.2 2.8])


% save the plot
% =============================================
plotEPS = sprintf('eps/%s.eps',dataFileName);
plotPNG = sprintf('png/%s.png',dataFileName);
print('-depsc','-r1200',plotEPS);
print('-dpng',plotPNG);
hold off;
